function data2fftFebSTHarmonic
format long;

febx50=load('febx50.out');
febx100=load('febx100.out');
febx150=load('febx150.out');
febx200=load('febx200.out');
febx250=load('febx250.out');
febx300=load('febx300.out');
febx350=load('febx350.out');
febx400=load('febx400.out');
febx450=load('febx450.out');
febx500=load('febx500.out');
febx550=load('febx550.out');
febx600=load('febx600.out');
febx650=load('febx650.out');

datax50=[febx50(:,1) febx50(:,2)];
datax100=[febx100(:,1) febx100(:,2)];
datax150=[febx150(:,1) febx150(:,2)];
datax200=[febx200(:,1) febx200(:,2)];
datax250=[febx250(:,1) febx250(:,2)];
datax300=[febx300(:,1) febx300(:,2)];
datax350=[febx350(:,1) febx350(:,2)];
datax400=[febx400(:,1) febx400(:,2)];
datax450=[febx450(:,1) febx450(:,2)];
datax500=[febx500(:,1) febx500(:,2)];
datax550=[febx550(:,1) febx550(:,2)];
datax600=[febx600(:,1) febx600(:,2)];
datax650=[febx650(:,1) febx650(:,2)];

NNN=length(febx50);

fftx50=fft(datax50(:,2))/NNN;
fftx100=fft(datax100(:,2))/NNN;
fftx150=fft(datax150(:,2))/NNN;
fftx200=fft(datax200(:,2))/NNN;
fftx250=fft(datax250(:,2))/NNN;
fftx300=fft(datax300(:,2))/NNN;
fftx350=fft(datax350(:,2))/NNN;
fftx400=fft(datax400(:,2))/NNN;
fftx450=fft(datax450(:,2))/NNN;
fftx500=fft(datax500(:,2))/NNN;
fftx550=fft(datax550(:,2))/NNN;
fftx600=fft(datax600(:,2))/NNN;
fftx650=fft(datax650(:,2))/NNN;

% Write the initial conditions to a file.
% These initial conditions are needed by the Fortran solvers.

ZeroModeLoc=201;
deltak=1/(datax50(NNN,1)-datax50(1,1));
kvec=transpose(0:deltak:(NNN-1)*deltak);
k0=1.52804650398968373
epsilon=2*k0*abs(fftx50(ZeroModeLoc))

% Compute M norms and exponential fit.

x50Mnorm=0;
x100Mnorm=0;
x150Mnorm=0;
x200Mnorm=0;
x250Mnorm=0;
x300Mnorm=0;
x350Mnorm=0;
x400Mnorm=0;
x450Mnorm=0;
x500Mnorm=0;
x550Mnorm=0;
x600Mnorm=0;
x650Mnorm=0;

for j=151:251
    x50Mnorm=x50Mnorm+abs(fftx50(j))^2;
    x100Mnorm=x100Mnorm+abs(fftx100(j))^2;
    x150Mnorm=x150Mnorm+abs(fftx150(j))^2;
    x200Mnorm=x200Mnorm+abs(fftx200(j))^2;
    x250Mnorm=x250Mnorm+abs(fftx250(j))^2;
    x300Mnorm=x300Mnorm+abs(fftx300(j))^2;
    x350Mnorm=x350Mnorm+abs(fftx350(j))^2;
    x400Mnorm=x400Mnorm+abs(fftx400(j))^2;
    x450Mnorm=x450Mnorm+abs(fftx450(j))^2;
    x500Mnorm=x500Mnorm+abs(fftx500(j))^2;
    x550Mnorm=x550Mnorm+abs(fftx550(j))^2;
    x600Mnorm=x600Mnorm+abs(fftx600(j))^2;
    x650Mnorm=x650Mnorm+abs(fftx650(j))^2;
end
 
positions=epsilon^2*k0*[0,50,100,150,200,250,300,350,400,450,500,550,600]';
xdata=epsilon^2*k0*linspace(0,650,100)';
Mnorms=(k0/epsilon)^2*[x50Mnorm,x100Mnorm,x150Mnorm,x200Mnorm,x250Mnorm,x300Mnorm,x350Mnorm,x400Mnorm,x450Mnorm,x500Mnorm,x550Mnorm,x600Mnorm,x650Mnorm]';
logMnorms=log(Mnorms);
coeffs=polyfit(positions,logMnorms,1);
Mpred=exp(polyval(coeffs,xdata));
delta=-coeffs(1)/2
L2Scaling=exp(2*delta.*(positions));
ScaledMnorms=Mnorms.*L2Scaling;

figure(3)
subplot(2,1,1)
hold off
scatter(positions,Mnorms,'filled','k')
hold on
plot(xdata,Mpred,'r','LineWidth',2)
text(0.05,0.15,strcat('Computed \delta=',num2str(delta)),'FontSize',13,'FontName','Times')
ax=gca;
ylabel('$\mathcal{M}$','Interpreter','Latex')
legend('data','best fit exp','Location','NorthEast')
title('$\mathcal{M}$','Interpreter','Latex')

subplot(2,1,2)
hold off
scatter(positions,ScaledMnorms,'filled','k')
hold on
plot([0 0.2],[(k0/epsilon)^2*x50Mnorm (k0/epsilon)^2*x50Mnorm],'r','LineWidth',2)
text(0.05,0.36,strcat('Using \delta=',num2str(delta)),'FontSize',13,'FontName','Times')
xlabel('$\chi$','Interpreter','Latex')
ylabel('$\mathcal{M}$','Interpreter','Latex')
leg1=legend('scaled data','$\mathcal{M}(0)$','Location','NorthEast');
set(leg1,'Interpreter','Latex');
title('Scaled $\mathcal{M}$','Interpreter','Latex')

